Práctica Parcial 3 - Teoría de la Información

Can Bacab Héctor Arturo, 150300161

Chi Polanco Adrián de Jesús, 150300103

Peralta Ramos José Francisco, 150300163

2020-12-04

El concepto de entropía

La entropía se puede interpretar como un indicativo de la complejidad en una serie de datos aleatorios \(\{X_1, X_2, \ldots, X_N\}\). La entropía usualmente se calcula utilizando una función de masa de probabilidad \(p_j\) que cumple con las siguientes propiedades:

  • \(p_j \ge 0, \, j \in 0,1, \ldots, N\)
  • \(\sum_{j=1}^N{p_j} = 1\)

R tiene múltiples paquetes y funcionalidades que permiten estimar la pmf de un conjunto de datos como el descrito anteriormente. El histograma es una herramienta que permite estimar la pmf de un conjunto de datos. En el ejemplo que sigue se muestra la forma de estimar la pmf en R:

set.seed(1234)                       # Para hacer el análisis reproducible
datos     <- rnorm(512,0,1)          # Se generan 512 valores normales
histogram <- hist(datos, plot=FALSE) # Se calcula el histograma
pmf       <- histogram$counts/sum(histogram$counts)  # Se calcula la pmf
sum(pmf)                             # Se verifica que cumpla con las propiedades
[1] 1

Ejercicios

  • Estimar la pmf utilizando el utilizando los paquetes ASH y KernSmoooth.
library(ash)
a <- ash1(bin1(datos, ab = c(min(datos), max(datos))))
[1] "ash estimate nonzero outside interval ab"
plot(a, type = "l")

sum(diff(a$x) * (head(a$y, -1) + tail(a$y, -1)))/2
[1] 0.9973948
library(KernSmooth)
k <- bkde(x = datos, range.x =  c(min(datos), max(datos)))
plot(k, type = 'l')

  • ¿Cuál es la ventaja de utilizar los métodos anteriores sobre el histograma?

El histograma es el método más viejo y menos sofisticado para estimar la densidad. ASH utiliza el método de average shifted histogram, el cual es más suave que el histograma y evita la sensibilidad a la elección del origen, pero sigue siendo computacionalmente eficiente.

KernSmooth utiliza el método de kernel density estimation. El kernel es una función simétrica, generalmente positiva, que se integra en 1. El enfoque de kernel density estimation supera la discreción de los enfoques del histograma al centrar una función de kernel suave en cada punto de datos y luego sumar para obtener una estimación de densidad.

  • Utilizando el comando hist y los paquetes ASH y KernSmooth verifique el tiempo requerido para estimar la densidad de una serie de datos Gaussianos con \(\mu=1\) y varianza \(\sigma^2=1\) y longitudes \(N=2^i, \, i=8,9,10, 11, \ldots 16.\) (es necesario incluir un gráfico en highcharter)
dataSeries <- list()
for (i in 8:16){
    dist <- rnorm(2^i, mean = 1, sd = 1)
    dataSeries[[(i-7)]] <- dist
}

dataAll = do.call(rbind, dataSeries)
tiempoHist <- list(NA, 9)
for (i in 1:9){
    start <- Sys.time()
    hist(dataAll[i,], plot=FALSE)
    end <- Sys.time()
    runtime <- end - start
    
    tiempoHist[[i]] <- runtime
}
th <- do.call(rbind, tiempoHist)
nombre <- rep("Hist", 9)
n_i <- (8:16)
th <- cbind(th, nombre, n_i)
tiempoASH <- list(NA, 9)
for (i in 1:9){
    start <- Sys.time()
    ash1(bin1(dataAll[i,], ab = c(min(dataAll[i,]), max(dataAll[i,]))))
    end <- Sys.time()
    runtime <- end - start
    
    tiempoASH[[i]] <- runtime
}
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
[1] "ash estimate nonzero outside interval ab"
tash <- do.call(rbind, tiempoASH)
nombre <- rep("ASH", 9)
n_i <- (8:16)
tash <- cbind(tash, nombre, n_i)
tiempoKS <- list(NA, 9)
for (i in 1:9){
    start <- Sys.time()
    bkde(x = dataAll[i,], range.x =  c(min(dataAll[i,]), max(dataAll[i,])))
    end <- Sys.time()
    runtime <- end - start
    
    tiempoKS[[i]] <- runtime
}
tks <- do.call(rbind, tiempoKS)
nombre <- rep("KS", 9)
n_i <- (8:16)
tks <- cbind(tks, nombre, n_i)
tiempos <- rbind(th, tash, tks)
tiempos.df <- as.data.frame(tiempos)
names(tiempos.df)[1] <- "time"
names(tiempos.df)[2] <- "method"
names(tiempos.df)[3] <- "i"
tiempos.df$time <- as.numeric(tiempos.df$time)
head(tiempos.df)
        time method  i
1 0.01270819   Hist  8
2 0.01890492   Hist  9
3 0.02113605   Hist 10
4 0.02066207   Hist 11
5 0.01725292   Hist 12
6 0.01745105   Hist 13
library(highcharter)

Gráfica de tiempos

hchart(tiempos.df, "column", hcaes(x = i, y = time, group = method))

La entropía utilizando el histograma

Volviendo de nuevo al ejemplo anterior, podemos estimar la entropía de Shannon, utilizando la pmf obtenida mediante el histograma y así obtener un estimador empírico de la entropía de Shannon. A continuación mostramos la forma de obtener la entropía de una serie de datos obtenida en ventanas independientes o contiguas de longitud \(512\):

set.seed(1234)
datos        <- rnorm(32768)
wLength      <- 512
plot(datos, type = "l", main="Serie de datos normal", ylab="Valores", xlab="Tiempo, t")

noVentanas    <- length(datos)/wLength
entropies     <- numeric(noVentanas)
index         <- numeric(noVentanas)
for(i in 1:noVentanas)
{
  
  dataW        <- datos[wLength*(i-1)+1:wLength*i]
  histo        <- hist(dataW, breaks=8,plot=FALSE)
  pmf          <- histo$counts/sum(histo$counts)
  entropies[i] <- -1*sum(pmf*log(pmf))
  index[i]     <- wLength*(i-1)+1
}
plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")

Preguntas

  1. ¿Porqué existen valores discontinuos en la entropía?

Porque en algunos de los resultados de la pmf el resultado es 0 y al calcular la entropía obtenemos valores NaN.

  1. ¿Con qué código soluciona el problema de las discontinuidades?

Podemos excluir los valores NaN con el parámetro na.rm

entropia <- function(datos, wLength)
    { 
    noVentanas    <- length(datos)/wLength
    entropies     <- numeric(noVentanas)
    index         <- numeric(noVentanas)
    for(i in 1:noVentanas)
    {
  
      dataW        <- datos[wLength*(i-1)+1:wLength*i]
      histo        <- hist(dataW, breaks=8,plot=FALSE)
      pmf          <- histo$counts/sum(histo$counts)
      entropies[i] <- -1*sum(pmf*log(pmf), na.rm = TRUE)
      index[i]     <- wLength*(i-1)+1
    }
    plot(index, entropies, type = "l", main="Entropías empíricas para datos normales", xlab="Tiempo, t", ylab="Valores de entropía")
    }
entropia(datos, wLength)

  1. Ahora calcule la entropía de una serie de datos normales, denotados por \(X_t\) (con \(\mu=0\) y \(\sigma=2\)) pero ahora añadanle (súmenle) una segunda función \(r_t\), es decir, hallen la entropía de la serie \(Y_t = X_t+r_t\) con \(r_t\) definida por: \[ r_t = \begin{cases} \sigma/4 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)

  1. Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma/2 & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)

  1. Repitan el paso \(3\) para \(r_t\) dada por: \[ r_t = \begin{cases} \sigma & t\ge 16384\\ 0 & \mbox{otro caso} \end{cases} \]
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512
entropia(datos, wLength)

  1. ¿Tiene algún efecto la longitud del salto en \(r_t\) en la forma de la entropía? Explique.

Entre mayor sea el valor de sigma se presentan menos picos en la entropía.

  1. ¿Qué sucede ahora si \(r_t\) es de la forma: \[ r_t = \begin{cases} \sigma & 15872 \le t\le 16896\\ 0 & \mbox{otro caso} \end{cases} \]?
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, 2)
wLength <- 512
entropia(datos, wLength)

Los picos de los valores intermedios disminuyen.

  1. Repita los pasos 3-7 pero ahora usando la entropía de Harvda con parámetro \(\alpha=3\) y \(\alpha=9\). ¿Qué efecto tiene \(\alpha\)?
havrda <- function(pmf, alpha) sum(pmf^alpha - 1) / (2^(1 - alpha) - 1)
entropiaHavrda <- function(datos, wLength, alphaVal)
    { 
    noVentanas    <- length(datos)/wLength
    entropies     <- numeric(noVentanas)
    index         <- numeric(noVentanas)
    for(i in 1:noVentanas)
    {
  
      dataW        <- datos[wLength*(i-1)+1:wLength*i]
      histo        <- hist(dataW, breaks=8,plot=FALSE)
      pmf          <- histo$counts/sum(histo$counts)
      entropies[i] <- havrda(pmf, alpha = alphaVal)
      index[i]     <- wLength*(i-1)+1
    }
    plot(index, entropies, type = "l", xlab="Tiempo, t", ylab="Valores de entropía")
    }
  • Paso 3
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512

\(\alpha = 3\)

entropiaHavrda(datos, wLength, 3)

\(\alpha = 9\)

entropiaHavrda(datos, wLength, 9)

  • Paso 4
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512

\(\alpha = 3\)

entropiaHavrda(datos, wLength, 3)

\(\alpha = 9\)

entropiaHavrda(datos, wLength, 9)

  • Paso 5
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512

\(\alpha = 3\)

entropiaHavrda(datos, wLength, 3)

\(\alpha = 9\)

entropiaHavrda(datos, wLength, 9)

  • Paso 7
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
wLength <- 512

\(\alpha = 3\)

entropiaHavrda(datos, wLength, 3)

\(\alpha = 9\)

entropiaHavrda(datos, wLength, 9)

Podemos observar que entre menor sea el valor de \(\alpha\) el valor de la entropía disminuirá.

Entropía utilizando ventanas deslizantes

El cálculo de la entropía por ventanas independientes dada arriba resulta útil en casos en dónde la función no tiene dependencia en los valores futuros (descorrelacionadas). Para el caso de funciones correlacionadas, el cálculo de la entropía por ventanas deslizantes traslapadas resulta útil para descubrir ciertas fenomenologías en los datos. El cálculo por ventanas deslizantes de tamaño \(W\) se realiza sobre una secuencia de datos \(X_1, X_2, \ldots, X_N\). La ventana (\(W\le N\)) se va deslizando sobre los datos con factor \(\Delta\) y de esta forma subconjuntos de los datos \(X_i\) toman la siguiente forma: \[ X(m; W, \Delta) = x_j \times \Pi(\frac{t-m\Delta}{W}-\frac{1}{2}), \] donde \(m\Delta \le j \le m\Delta + W\) y \(m=0,1,2, \ldots\). Finalmente se puede graficar \(nW + \Delta, n=1,2,3, \ldots\) contra las entropías y verificar algún patrón en los datos.

Ejercicios

  • Implementar en R la metodología del cálculo de la entropía de Harvda normalizada por ventanas deslizantes. La función debe tener la forma harvda_deslizante(datos, w.length=512, s.factor=10, a.parameter=0.8, ent.type=c("hist", "ash", "kern")), donde w.length es la longitud de la ventana, s.factor es el factor de deslizamiento y a.parameter es el parámetro \(\alpha\) de la entropía de Harvda. Además la función puede calcular la entropia usando el histograma, por el método ash o por alguna metodología kernel (con el parámetro ent.type).

  • Aplicar la entropía calculada con los datos generados anteriormente, es decir:

  • Los pasos \(3,4,5\) y \(7\) en donde los datos se dan como \(X_t+r_t\).

havrdaNormalizada <- function(pmf, N, a.parameter) 
    {
    sum(pmf^a.parameter - 1) / (N^(1 - a.parameter) - 1)
    }
havrdaDeslizante <- function(datos, w.length = 512, s.factor = 10, a.parameter = 0.8, ent.type=c("hist", "ash", "kern"))
    {
    N <- length(datos)
    m <- 0
    n <- 1
    
    inf <- 0
    sup <- 0
    
    entropies <- c()
    index <- c()
    
    inf <- m * s.factor + 1
    sup <- m * s.factor + w.length
    
    while (sup <= N)
        {
        dataW <- datos[inf:sup]
        pmf <- switch (ent.type,
            "hist" = {hist <- hist(dataW, breaks = 1000, plot = FALSE)
                      hist$counts / sum(hist$counts)},
            "ash" = {invisible(capture.output(hist <- ash1(bin1(dataW))))
                     hist$y / sum(hist$y)},
            "kern" = {hist <- bkde(dataW)
                      hist$y / sum(hist$y)})
        
        entropies <- append(entropies, havrdaNormalizada(pmf, N, a.parameter))
        index <- append(index, n*w.length + s.factor)
        
        m <- m + 1
        n <- n + 1
        
        inf <- m * s.factor + 1
        sup <- m * s.factor + w.length
        }
    
    return(list(index, entropies))
    }
  • Paso 3
library(zeallot)
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/4, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist") 
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash") 
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

  • Paso 4
r_t <- function(t, sigma) ifelse(t >= 16384, sigma/2, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist") 
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash") 
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

  • Paso 5
r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist") 
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash") 
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

  • Paso 7
r_t <- function(t, sigma) ifelse(15872 <= t & t<= 16896, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)
c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "hist") 
plot(index, entropies, type = "l", main="Usando hist", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "ash") 
plot(index, entropies, type = "l", main="Usando ash", xlab="Tiempo, t", ylab="entropy")

c(index, entropies) %<-% havrdaDeslizante(datos, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

  • ¿Tiene algún efecto el parámetro \(a\) en la forma de la entropía? ¿Tiene alguna ventaja calcular la entropía de Harvda por otro método diferente al histograma?

Comparamos utilizando kern y modificando el valor de \(\alpha\)

r_t <- function(t, sigma) ifelse(t >= 16384, sigma, 0)
set.seed(1234)
datos <- rnorm(32768, mean = 0, sd = 2) + r_t(1:32768, sigma = 2)

\(\alpha\) = 1.5

c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 1.5, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

\(\alpha\) = 3

c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 3, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

\(\alpha\) = 5

c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 5, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

\(\alpha\) = 7

c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 7, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

\(\alpha\) = 9

c(index, entropies) %<-% havrdaDeslizante(datos, a.parameter = 9, ent.type = "kern") 
plot(index, entropies, type = "l", main="Usando kern", xlab="Tiempo, t", ylab="entropy")

\(\alpha\) traslada la función sobre el eje y también disminuye el tamaño de la función.

  • Además del histograma y los estimadores tipo kernel, existen otros métodos para estimar la distribución de una serie de datos. Investigue: ¿en qué consiste la entropía de permutación?

Según Judge y Henry (2019) la entropía de permutación es una herramienta robusta de series de tiempo que proporciona una medida de cuantificación de la complejidad de un sistema dinámico al capturar las relaciones de orden entre los valores de una serie de tiempo y extraer una distribución de probabilidad de los patrones ordinales.

Referencias

Deng, H. & Wickham, H.. (2011). Density estimation in R.

Henry, M. (2020). Permutation Entropy. Aptech

Barnier, J. (2020). Rmdformats: HTML Output Formats and Templates for Rmarkdown Documents. https://github.com/juba/rmdformats.